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Abstract 

Let fi,...,f a £ K[xi, . . . , x m ] be a system of polynomials generating a zero-dimensional 
ideal I, where K is an arbitrary algebraically closed field. Assume that the factor algebra 
A — M.[xi, . . . ,x m ]/I is Gorenstein and that we have a bound 8 > such that a basis for A 
can be computed from multiples of fi, . . . , f 3 of degrees at most 8. We propose a method using 
Sylvester or Macaulay type resultant matrices of fi, . . . , f s and J, where J is a polynomial of 
degree 8 generalizing the Jacobian, to compute moment matrices, and in particular matrices of 
traces for A. These matrices of traces in turn allow us to compute a system of multiplication 
matrices {M^Ji = 1, . . . , ra} of the radical \ZX, following the approach in the previous work by 
Janovitz-Freireich, Ronyai and Szanto. Additionally, we give bounds for 8 for the case when I 
has finitely many projective roots in PJj . 

Keywords: Moment Matrices; Matrices of Traces; Radical Ideal; Solving polynomial systems 

1 Introduction 

This paper is a continuation of our previous investigation in [22\ 123] to compute the approximate 
radical of a zero dimensional ideal which has zero clusters. The computation of the radical of a zero 
dimensional ideal is a very important problem in computer algebra since a lot of the algorithms for 
solving polynomial systems with finitely many solutions need to start with a radical ideal. This is 
also the case in many numerical approaches, where Newton- like methods are used. From a symbolic- 
numeric perspective, when we are dealing with approximate polynomials, the zero-clusters create great 
numerical instability, which can be eliminated by computing the approximate radical. 

The theoretical basis of the symbolic-numeric algorithm presented in j22, : 23] was Dickson's lemma 
|14j . which, in the exact case, reduces the problem of computing the radical of a zero dimensional 
ideal to the computation of the nullspace of the so called matrices of traces (see Definition fT4|): in 
|22[I23| we studied numerical properties of the matrix of traces when the roots are not multiple roots, 
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but form small clusters. Among other things we showed that the direct computation of the matrix of 
traces (without the computation of the multiplication matrices) is preferable since the matrix of traces 
is continuous with respect to root perturbations around multiplicities while multiplication matrices 
are generally not. 

It turns out that the computationally most expensive part of the method in |22l I23| is the com- 
putation of the matrix of traces. We address this problem in the present paper, and give a simple 
algorithm using only Sylvester or Macaulay type resultant matrices and elementary linear algebra to 
compute matrices of traces of zero dimensional ideals satisfying certain conditions. 

More precisely, we need the following assumptions: let f = [fi, ■ ■ ■ , f s ] be a system of polynomials 
of degrees d\ > • • • > d s in K[x], with x = [xi, . . . , x rn ], generating an ideal T in K[x], where K is an 
arbitrary algebraically closed field. We assume that the algebra A := K[x]/Z is finite dimensional over 
K and that we have a bound 8 > such that a basis S = [bi , . . . , 6jv] of A can be obtained by taking 
a linear basis of the space of polynomials of degree at most S factored by the subspace generated by 
the multiples of /i, . . . , f s of degrees at most 6. By slight abuse of notation we denote the elements of 
the basis S which are in A and some fixed preimages of them in K[x] both by b\, . . . , &j\r. Thus we can 
assume that the basis S consists of monomials of degrees at most 5. Note that we can prove bounds 
6 = Y^X 1 di ~ 171 ( or ^ = Eti di — m if s — m) if I has only finitely many projective common roots 
in Pjp and have no common roots at infinity, using a result of Lazard [33] (see Theorem [3]). 

Furthermore, we also assume that A is Gorenstein over K (see Definition [I} . Note that in practice 
we can easily detect if A is not Gorenstein (see Remark ITU]) . Also, a random change of projective 
variables can eliminate roots at infinity with high probability when they are in finite number, but we 
will address the necessity of this assumption in an upcoming paper. 

The main ingredient of our method is a Macaulay type resultant matrix MacA(f ), which is defined 
to be a maximal row-independent submatrix of the transpose matrix of the degree A Sylvester map 
(<ji, . . . ,g s ) i — > J2i=i hdi <= K[x]a for A < 28 + 1 (see Definition EJ. Using our assumptions on A, 
we can compute a basis S of A using MacA(f), and we also prove that a random element y of the 
nullspace of MacA(f) provides a non-singular N x N moment matrix 9)1,5 (y) with high probability 
(similarly as in [31jV This moment matrix allows us to compute the other main ingredient of our 
algorithm, a polynomial J of degree at most <5, such that J is the generalization of the Jacobian of 
/i, . . . f s in the case when s = in. The main result of the paper now can be formulated as follows: 

Theorem Let S = [b\, . . . , i>jv] be a basis of A with deg(6j) < 5. With J as above, let Syl iS (J) be 
the transpose matrix of the map $3 i=1 Cibi *— > J • Xa=i K[x]a for Ci € IK. Then 

[rr(6i6i)C =1 = Syl s (J)-X, 

where X is the unique extension of the matrix 9Jts(y) such that MacA(f) ■ X = 0. 

Once we compute the matrix of traces R := [Tr(bibj)]fj =1 and the matrices R Xk := [Tr(xkbibj)]^j =1 
Syl s (xfc J) -X for k = 1, . . . , m, we can use the results of [2HI23] to compute a system of multiplication 
matrices for the (approximate) radical of X as follows: if R is a (numerical) maximal non-singular 
submatrix of R and R Xk is the submatrix of R Xk with the same row and column indices as in R, then 
the solution M Xk of the linear matrix equation 

RM Xk = R Xk 
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is an (approximate) multiplication matrix of Xk for the (approximate) radical of X. See |23j for the 
definition of (approximate) multiplication matrices. Note that a generating set for the radical VX can 
be obtained directly from the definition of multiplication matrices, in particular, it corresponds to the 
rows of the matrices M Xl , . . . , M Xm . 

We also point out that in the s = m case these multiplication matrices M Xk can be obtained even 
more simply using the nullspacc of MacA(f) and the Jacobian J of f , without computing the matrices 
of traces. 

We also note here that in a follow up paper we will consider an extension of our present results 
which works also in the non-Gorenstein case to compute the matrices of traces. Furthermore, that 
paper will also extend our results to the affine complete intersection case using Bczout matrices. 

2 Related Work 

The motivation for this work was the papers [311 E2] where they use moment matrices to compute 
the radical of real and complex ideals. They present two versions of the method for the complex case: 
first, in [32] they double up the machinery for the real case to obtain the radical of the complex ideal. 
However, in |31j they significantly simplify their method and show how to use moment matrices of 
maximal rank to compute the multiplication matrices of an ideal between X and its radical s/T. In 
particular, in the Gorenstein case they can compute the multiplication matrices of X. In fact, in [31] 
they cite our previous work [22| to compute the multiplication matrices of VT from the multiplication 
matrices of X, but the method proposed in the present paper is much simpler and more direct. 

Note that one can also obtain the multiplication matrices of X with respect to the basis S = 
[&i,...,6/v] by simply eliminating the terms not in S from x^bi using Mac<s+i(f). The advantage 
of computing multiplication matrices of the radical y/T is that it returns matrices which are always 
simultaneously diagonalizablc, and possibly smaller than the multiplication matrices of X, hence easier 
to work with. Moreover, if S contains the monomials 1, x\, . . . , x m , one eigenvector computation yield 
directly the coordinates of the roots. 

Computation of the radical of zero dimensional complex ideals is very well studied in the literature: 
methods most related to ours include [HE] where matrices of traces are used in order to find generators 
of the radical, and the matrices of traces are computed using Grobner Bases; also, in [1] they use the 
traces to give a bound for the degree of the generators of the radical and use linear solving methods 
from there; in [19j they describe the computation of the radical using symmetric functions which 
are related to traces. One of the most commonly quoted method to compute radicals is to compute 
the projections X n K[xj] for each i = 1, . . . , m and then use univariate squarefree factorization (see 
for example [171 [26] [TO] [20] ). The advantage of the latter is that it can be generalized for higher 
dimensional ideals (see for example [25]). We note here that an advantage of the method using matrices 
of traces is that it behaves stably under perturbation of the roots of the input system, as was proved in 
[2"3] . Other methods to compute the radical of zero dimensional ideals include [531 EE HE1 12H1 1301 132] • 
Applications of computing the radical include |21j , where they show how to compute the multiplicity 
structure of the roots of X once the radical is computed. 

Methods for computing the matrix of traces directly from the generating polynomials of X, without 
using multiplication matrices, include [H [6] where they use Newton Sums, [3 [9] where they use 
residues and [12] using resultants. Besides computing the radical of an ideal, matrices of traces have 
numerous applications mainly in real algebraic geometry [2j 1351 [4] , or in |36j where trace matrices are 
applied to find separating linear forms detcrministically. 
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3 Moment Matrices and Matrices of Traces 

Let f = [/i, . . . , f a ] be a system of polynomials of degrees d\ > ■ • • > d s in K[x], where x = [x%, . . . , x m ] 
and IK is an arbitrary algebraically closed field. Let X be the ideal generated by /i, . . . , f s in K[x] and 
define A := K[x]/X. We assume throughout the paper that A is a finite dimensional vector space over 
K and let A* denote the dual space of A. 

Let us first recall the definition of a Gorenstein algebra (c.f. [UJ 133 [T5J [31]). Note that these 
algebras are also referred to as Frobenius in the literature, see for example [2j. 

Definition 1 A finite dimensional K- algebra A is Gorenstein (overWL) if there exists a nondegenerate 
'K-bilinear form B(x,y) on A such that 

B(ab, c) = -B(a, be) for every a, b, c G A. 

Note that this is equivalent to the fact that A and A* are isomorphic as A modules. It is also 
equivalent to the existence of a K-lincar function A : A — > K such that the bilinear form B(a, b) := 
A(ab) is nondegenerate on A. 

Assumption 2 Throughout the paper we assume that A is Gorenstein. Furthermore, we also assume 
that we have a bound S > such that 

N := dim K K[x] 5 /(/ 1 , ...,/.}* = dim K K[x] d /(/ 1 , . . . , f s ) d (1) 

for all d > S and that 

N = dimA (2) 

Here 

ii /-..<: |X f*Pi : de S^) <d~d^. (3) 

We fix S = [&i, . . . , &at] a monomial basis for A such that deg(bi) < 5 for all i = 1, ... ,7V. Lei D be 
the maximum degree of the monomials in S . Thus D < 5. 

We have the following theorem giving bounds for 5 in the case when f has finitely many projective 
roots. 

Theorem 3 Let f = [fi, . . . , f s ] be a system of polynomials of degrees d\ > ■ ■ ■ > d s in K[x]. Assume 
that fx, . . . , f s has finitely many projective common roots in PJj? . Assume further that fi, f%, ■ ■ ■ , f s 
have no common roots at infinity. Then: 

1. If s — m then for S := X)"=i(^» — 1) conditions (QJ) and (0) are satisfied. Furthermore, in this 
case A is always Gorenstein. 

2. If s > m then for 5 := Y^i^ 1 di — m conditions (Qp and (0) are satisfied. 
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Proof. 

For the first assertion let f h be the homogenization of f using a new variable x m +\. Using our 
assumption that f h has finitely many roots in PJ^ and s = m, one can see that (f h ) is a regular 
sequence in R := K[i'i, . . . , x m , x m +i}- 

Define the graded ring B := R/(i h ). Following the approach and notation in [35|, we can now 
calculate the Hilbcrt series of B, defined by H(B, A) = J^d ^B{d)X d , where Hb is the Hilbcrt function 
of B. We have 

H(R X) — H{B > X) 

and using the simple fact that 

H(R, A) = — — - 

v 1 (1 - A) m+1 

we obtain that 

(l + A + .-. + A^-^-.q + A + .-. + A^- 1 ) 
H(B, A) = ^— 

= g(X)(l + A + ...), 

where 

g(X) = (1 + A + • • • + A*- 1 ) • • • (1 + A + • • • + A dm_1 ). 
This implies that the Hilbcrt function 

H B (8) = H B (S + 1)=Hb{5 + 2) = ... 

Note that dehomogenization induces a linear isomorphism Bd — > K[x]d/(/i, • • • , f s )d, where Bd stands 
for the degree d homogeneous part of B. From this, using that there are no common roots at infinity, 
we infer that for d > 6 dimuK^xJci/ (fi, . . . , f s )d = dims .4 = N, which implies |T]) and ([2|). 
Note that the common value N = TLb(5) is the sum of the coefficients of g, which is 

i=l 

To prove that A is Gorenstein, we cite [151 Proposition 8.25, p. 221] where it is proved that if 
/ij---j/m is an affinc complete intersection then the Bezoutian Bij 1 _,_j m defines an isomorphism 
between A* and A. 

To prove the second assertion we note that [331 Theorem 3.3] implies that 

dim K Bs = dim K B s+1 = 

From here we obtain (Q} and ^ as in the Case 1. 

□ 



Remark 4 JVote that in general Id ^ (fi, ■ ■ ■ , f s )d, where Id is the set of elements of I with degree 
at most d and (/i, . . . ,f s )d was defined in |3p. This can happen when the system has a root at 
inhnity, for example, if fx = x + 1, f2=x then Jq = span K (l) but (/i, /2)o = {0} However, using the 
homogenization fx, ■ ■ ■ ,fs, the degree d part of the homogenized ideal is always equal to the space 
spanned by the multiples of /{',..., f£ of degree d. The above example also demonstrates that dim„4 
is not always the same as dimK[x],j/(/i, . . . ,f s )d for large enough d, because above dim .4 = but 
dimK[a;, y] d /(fx, f 2 )d = 1 for all d > 0.' 
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Next we will define Sylvester and Macaulay type resultant matrices for fx , . . . f s . 
Definition 5 Define 

A := max(<5, 2D + 1) 

where S and D are defined in Assumption^ 

Let Syl A (f) be the transpose matrix of the linear map 

0K[x] A _ di — » K[x] A (4) 

i 

s 

(gi,...,g a ) i-> ^2,fi9i 

i=l 

written in the monomial bases. So, in our notation, Syl A (f) will have rows which correspond to all 
polynomials fiX a of degree at most A. 

Let MacA(f) be a row submatrix of Syl A (f) of maximal size with linearly independent rows. 

Remark 6 In the case where s = m, for generic f we can directly construct MacA(f) by taking the 
restriction of the map HP to 

m 

05, (A) -+K[x] A 

i=l 

where £;(A) = span{x" : |a| < A — d,, Vj < i, otj < dj}. 

Here MacA (f ) is a submatrix of the classical Macaulay matrix of the homogcnization of f and some 
fm+i , where f^ + \ is any homogeneous polynomial of degree A— S: we only take the rows corresponding 
to the polynomials in f. Since the Macaulay matrix is generically non-singular, MacA(f) will also be 
generically full rank. 

Note that with our assumption that fi , . . . , f m has hnitely many projective roots, we have that 
MacA(f) has column corank TV := YliLi di- 

Since A > 5, by Assumption [5] the corank of MacA(f) = TV, where TV is the dimension of A. Also, 
we can assume that the elements of the basis S of A are monomials of degree at most 5, and that the 
first columns of MacA(f) correspond to the basis S of A. 

Fix an element 

y = [y a : a e N m , \a\ < A] T 
of the nullspacc Null(MacA(f)), i-e. Mac A (f) ■ y = 0. 

Definition 7 Let S be the basis of A as above, consisting of monomials of degree at most D. Using 
y we can define A y € A* by A y (g) := ^2 x a e gy a ga, where g — X)x°es-9 QxQ e A- Note that every 
A G A* can be defined as A y for some y € Null(Mac A (f)) or more generally with an element o/K[x]* 
which vanishes on the ideal T. 

Define the moment matrix 9Jts(y) to be the TV x TV matrix given by 

where a and (3 run through the exponents of the monomials in S . Note that SUI5 is only a submatrix 
of the usual notion of moment matrices in the literature, see for example Ml] . 

For p € A, we define the linear function p-A6 A* as p ■ A(g) := A(j>g) for all g 6 A. 
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Remark 8 If one considers a linear function A on A, such that the bilinear form (x, y) i— ► A(xy) is 
nondegenerate on A, then the moment matrix corresponding to this A will be the one whose (i,j)-th 
entry is just A (6, fry). Moreover, for g,h £ A 

A y (gh) = cocffs(. g ) T • Tt s (y) ■ coeS s (h) 

where coeffg(p) denotes the vector of coefficients ofpGA in the basis S. 

The following proposition is a simple corollary of [3TJ Prop 3.3 and Cor. 3.1]. 

Proposition 9 Let y be a random element of the vector space Null(MacA(f))- With high probability, 
yftsiy) is non-singular. 

Remark 10 Using the above proposition, one can detect whether the algebra A is not Gorenstein 
with high probability by simply computing the rank of 50? s (y) for (perhaps several ) random elements 
y in Null(Mac A (f))' 

Remark 11 By pffl Theorem 2.6 and Lemma 3.2] one can extend y toy £ K N such that the inhnitc 
moment matrix 9Jt(y) := [y a +0]a,/3eN m has the same rank as SDTs(y) and the columns of £Dt(y) vanish 
on all the elements of the ideal X. 

Next we define a basis dual to S = \b\, . . . , 6 at] with respect to the moment matrix 97ls(y). Using 
this dual basis we also define a polynomial J which is in some sense a generalization of the Jacobian 
of a well-constrained polynomial system. 

Definition 12 From now on we fix y € Null(MacA(f)) such that 9Jls(y) is invertible and we will 
denote by A the corresponding element A y £ A* . We define 

ms\y) =■■ M5=i- 

Let b* := YljLi c ji^j- Then [b^, . . . , b* N ] corresponds to the columns of the inverse matrix 9Jl^ 1 (y) 
and they also form a basis for A. Note that we have A(fe;6*) = 1, if i = j, and otherwise. 
Define the generalized Jacobian by 

N 

J := bib* mod X (5) 

i=l 

expressed in the basis S — [b±, . . ., 6jy] of A. 

Remark 13 Note that since Xw=i ^* has degree at most 2D, and A > 2D, we can use MacA(f) to 
End its reduced form, which is J. Because of this reduction, we have that deg( J) < D < 5. 

Also note that the notion of generalized Jacobian was also introduced in {2f. Its name come from 
the fact that if s = m and if A is the so called residue (c.f \15]). then bib* = J is the Jacobian 

of fi, ■ ■ ■ 

We now recall the definition of the multiplication matrices and the matrix of traces as presented 
in [23]. 
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Definition 14 Let p £ A. The multiplication matrix M p is the transpose of the matrix of the multi- 
plication map 

M p :A — >A 

g^pg 

written in the basis S. 

The matrix of traces is the N x N symmetric matrix: 

R=[Tr(b l b 3 )]f j=1 

where Tr(pq) := Tr{M pq ), M pq is the multiplication matrix of pq as an element in A in terms of the 
basis S = [61,..., 6jv] and Tr indicates the trace of a matrix. 

The next results relate the multiplication by J matrix to the matrix of traces R. 

Proposition 15 Let Mj be the multiplication matrix of J with respect to the basis S . We then have 
that 

Mj = [Tr(hb*)}» j=1 . 
Proof. Let A £ A* be as in Definition [T^l For any h £ A we have that 

N N 

h = J2A(hb j )b* = J2A(hb*)b j 

J = l 3=1 
N 



hbi = J2 A (hb*bi)b j = M h [i,j] = A(hb*bi 

N N 

Tr(h) = J2 Khb*bi) = A(h ]T &*&*). 



; = 1 



Since J = J^Zi in A wc have Tr(h) = A(hJ). 
Therefore 

Mj[i,j] = A(Jb*h) = Tr{b*h) = Tr(bib*) 



□ 



Corollary 16 

or equivalently J ■ A = Tr in A* . 



Mj ■ m s (y) = [Trihbj)}" = R, 



Proof. The coefficients of b* in the basis S = [bi, . . . , bpf] are the columns of 9Jl s 1 (y), which implies 
that 

Mj = [Tr&b*)]"^ = [Trihb^^-Ttfiy). 
Therefore we have that Mj ■ m s (y) = [Tripibj)]? , =1 . □ 



Finally, we prove that the matrix of traces R can be computed directly from the Sylvester matrix 
of f\, . . . , f s and J, without using the multiplication matrix Mj. First we need a lemma. 
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Lemma 17 There exists a unique matrix 5Ks(y) of size |Mon<(A) — S\ X \S\ such that 



Mac A (f) 



«s(y) 



Proof. By our assumption that the first columns of MacA(f) correspond to S we have 



Mac A (f) 



B 



A 



where the columns of B are indexed by the monomials in S. Note here that by Assumption [5] the 
rows of MacA(f) span 2a, and the monomials in S span the factor space K[x]a/2a- These together 
imply that the (square) submatrix A is invertible. 
Then 



B 



A 



-A^B 



which implies that 



Mac A (f) 



5Hs(y) 



where «R s (y) = -A~ l B • 9Jt s (y). 



□ 



By construction, the column of £0Ts(y) indexed by bj € S corresponds to the values of bj ■ A G ^4* 
on b\, . . . , 6jv- The same column in £Ks(y) corresponds to the values of bj ■ A on the complementary 
set of monomials of Mon<(A). The column in the stacked matrix corresponds to the value of bj ■ A 
on all the monomials in Mon<(A). To evaluate bj ■ A(p) for a polynomial p of degree < A, we simply 
compute the inner product of the coefficient vector of p with this column. 

Definition 18 Let S = [b\, . . . , bpj] be the basis of A as above, and let P G K[x] be a polynomial of 
degree at most D + 1 . 

Define Sylg(P) to be the matrix with rows corresponding to the coefficients of the polynomials 
(biP), . . . , (b^P) in the monomial basis Mon<(A) (we use here that deg(fr;) < D, thus deg(6jP) < 

2D + 1 < a;. 

Furthermore, we assume that the monomials corresponding to the columns of Sy\ s (P) are in the 
same order as the monomials corresponding to the columns of MacA(f). 
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Theorem 19 



Syl s (J) 



OT s (y) 



= [Tr(6A')£ = i 



Proof. Since the j-th column of the matrix 



3Ks(y) 



«s(y) 



represents the values of 6j • A on all the monomials of degree less than or equal to A, and the i-th 
row of Sylg(J) is the coefficient matrix of bi J, we have 



Syl s (J) 



ans(y) 



*s(y) 



= [{^ ■ A)(biJ)]^ j=1 



= [TrQnbi)]^. 



□ 



We can now describe the algorithm to compute a set of multiplication matrices M Xi , i = 1, . . . ,m 
of the radical yX of X with respect to a basis of 'K.\x\/\/T. To prove that the algorithm below is 
correct we need the following result from [22 Proposition 8.3] which is the consequence of the fact 
that the kernel of the matrix of traces corresponds to the radical of A: 

Proposition 20 Let R be a maximal non-singular submatrix of the matrix of traces R. Let r be the 
rank of R, and T := . . . ,b ir ] be the monomials corresponding to the columns of R. Then T is a 
basis of the algebra 'K[x]/y/T and for each k = 1, . . . , m, the solution M Xk of the linear matrix equation 

RM Xh = R Xk 

is the multiplication matrix of Xk for \fl with respect to T . Here R Xk is the r x r submatrix of 
[Tr(xkbibj)]^j =1 with the same row and column indices as in R. 



Algorithm 21 Input: f = • • • , / s ] € K[x] of degrees di, . . . ,d s generating an ideal T and 6 > 
such that they satisfy the conditions in Assumption^ An optional input is D < S, which by default 
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set to be 8. 



Output: A basis T for the factor algebra K[x]/vT and a set of multiplication matrices {M Xi \i 
I m} of \jX with respect to the basis T . 

1. Compute Mac A (f) for A := max(2D + 1, 6) 

2. Compute a basis S o/K[x]A/(f)A such that the polynomials in S have degrees at most D. Let 
S = [h, . . . ,b N ]. 

3. Compute a random combination y of the elements of a basis of NuII(MelCa({)). 

4- Compute the moment matrix SDts(y) defined in Definition^ and $Hs(y) defined in Lemma [77[ 

5. Compute ^^^^(y) and the basis [&£,..., b^] defined in Definition ]!^ 

6. Compute J = YliLi ^iK mod X using MacA(f)- 

7. Compute Sy\ s (J) and Sy\ s (xkJ) for k = 1, . . . ,m defined in Definition \18\ 

8. Compute 



R = [Tr(&!&i)]&=i 



and 



^ t :=trr(a! fc Mj)]5=i = 



Syl s (J) 



OT s (y) 



«s(y) 



Syl s (x k J) 



aJts(y) 



SRs(y) 



for k = 1, . . . ,m. 



9. Compute R, a maximal non-singular submatrix of R. Let r be the rank of R, and T := 
[bi 1 , . . . , bi r ] be the monomials corresponding to the columns of R. 

10. For each k = 1, . . . ,m solve the linear matrix equation RM Xk = R Xk , where R Xk is the submatrix 
of R Xk with the same row and column indices as in R. 

Remark 22 Since the bound given in Theorem might be too high, it seems reasonable to design 
the algorithm in an iterative fashion, similarly to the algorithms in f3H 1321 140f . in order to avoid 
nullspacc computations for large matrices. The bottleneck of our algorithm is doing computations 
with MacA(f), since its size exponentially increases as A increases. 

Remark 23 JVote that if s = m then we can use the conventional Jacobian of f\, . . . , f m in the place 
of J, and any |Mon<(A)| x \S\ matrix X such that it has full rank and MacA(f) ■ X = in the place 
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of 



«s(y) 



Even though this way we will not get matrices of traces, a system of multiplication matrices of the 
radical VT can still be recovered: if Q denotes a maximal non-singular submatrix of Sy\ s (J) ■ X , 
and Q Xk is the submatrix of Sy\ s (xk J) • X with the same row and column indices as in Q, then the 
solution M Xk of the linear matrix equation QM Xk = Q Xk gives the same multiplication matrix of \[X 
w.r.t. the same basis T as the above Algorithm. 

Remark 24 As M Xk is the matrix of multiplication by xu modulo the radical ideal \/T, its eigen- 
vectors are (up to a non-zero scalar) the interpolation polynomials at the roots of X. Similarly the 
eigenvectors of the transposed matrix M Xk are (up to a non-zero scalar) the evaluation at the roots 
(ofl (see f34l 115} for more details). The vector which represents this evaluation at £ in the dual 
space A* is the vector of values of [&i, ... , 6 at] at Q. To obtain these vectors, we solve the generalized 
eigenvalue problem (R Xk — zR 1 ) w = and compute v = R* w. The vectors v will be of the form 
[61(C) j ■ ■ ■ 1 frjv(C)] f° r C a root of J. If b\ = 1, 62 = Xi, ... , b m+ i = x rn , we can read directly the 
coordinates of £ from this vector. 



4 Examples 

In this section we present three examples. Each of them has three polynomials in two variables. The 
first one is a system which has roots with multiplicities, the second one is a system which has clusters 
of roots, and the third one is a system obtained by perturbing the coefficients of the first one. For each 
of them we compute the Macaulay matrix MacA(f ), the vector y in its nullspace, the moment matrix 
VJlg(y), the polynomial J, the matrix of traces R and the (approximate) multiplication matrices of 
the (approximate) radical, following Algorithm 1211 
The exact system: 



3x\ + 18x13:2 - 48iri + 21x2 ~ 114x2 + 156 



— 5X2 + 6x2 -f- X2 + 5 
|_ 81 2 _ 163 2 , 2: 
h "4"X 1 X2 4^ x l + ^ 

+4x2 + 2x2 + 3 



-xix 2 + 



XiX 2 + ^XiX 2 



f has common roots (—1,3) of multiplicity 3 and (2, 2) of multiplicity 2. 
The system with clusters: 



3xi + 17.4xiX2 - 46.5xi + 23.855x2 - 127.977x 2 + 171.933 
xl - 72.943x^X2 + 139.617x^ - 8.417xix^ - 124.161xix 2 

+295.0283xi - 5x^ + 6x^ + x 2 + 5 
x'l + 21.853x^X2 - 43.658x5 - 27.011xix| + 185.548xix 2 

-274.649x1 - xl + 4x^ + 2x 2 + 3 



f has two clusters: (-1,3), (-0.9,3), (-1.01,3.1) and (2,2), (1.9,2) each of radius 10" 
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The perturbed system: 



3x( + 18x12:2 - 48a;i + 21.001a:2 - 113.99922 + 156.001 
1. OOl^i - 64.751S1X2 + 123.250^1 - 152.750a;ia:2 

+605.751ajix 2 - 587.500a:i - 4. 999^2 + 6.0001:1:1 + x 2 + 5 
x\ + 20.249x^2 - 40.750^ + 5.249X1XJ + 21. 749x1X2 

-75.5xi - 1. 001x2 + 4x? + 2x 2 + 3 



is obtained from f by a random perturbation of size 1CP 3 . This system has no common roots. 

We set 8 = 6, D = 2 and A = 6. The Sylvester matrices in all three cases were size 28 x 28 
and in the first two cases they had rank 23 while in the last case it was full rank. In the first two 
cases the fact that the corank is 5 indicates that there are 5 solutions, counting multiplicities. For 
these cases we computed a basis 5 := [l,Xi,X2,XiX2,x{\ for the factor algebra by taking maximum 
rank submatrices of the Macaulay matrices. In the third case, we simply erased the columns of the 
Macaulay matrix corresponding to the monomials in S. From here, we chose random elements in the 
nullspaces of the (cropped) Macaulay matrices to compute the moment matrices: 



— 6 
- 7 34 



-6 
7 

II) 

Jo 

-36 



-he 



1.787 
-20.059 






-7.207 
1.219 
1.787 



O 

1.787 
1.219 
7.702 
-43.499 





-20.059 

1.787 
-43.499 
-43.644 



-7.428 

-0.858 1.428 
-4.848 -0.858 





-0.858 

1.428 
-5.719 
-5.125 





-4.848 
-0.858 
-5.125 
-39.404 



The polynomials J, computed from the moment matrices are: 

T K 3 26 1 1 2 

J = 5 "T0 :Cl "T5 X2 -30 XlX2 -5 a:i 
J = 5 + 0.916:ci - 1.952^2 - 0.636^i^ 2 - 0.106a: 2 

J = 4.999 - 0.306xi - 1.733x2 - 0.030x1X2 - 0.200x 2 . 



After computing the matrices Syl 5 ( J) and £Rs (y), we obtain the matrices of traces: 



4.999 
0.990 
13.100 
-1.031 
10.440 



0.990 
10.440 
-1.031 
23.812 
12.100 



5 0.995 

0.995 10.999 

13.002 -1.017 
-1.017 25.0256 

11.003 12.870 



13.100 
-1.031 
35.610 
-11.206 
23.812 

13.002 
-1.015 
35.013 
-11.061 
25.0519 



-1.031 10.440 

23.812 12.100 

-11.206 23.812 

56.533 21.337 

21.337 31.729 

-1.017 11.003 

25.019 12.913 

-11.064 25.029 

59.129 22.770 

22.644 34.968 



The first matrix R has rank 2, while R and R have rank 5. In the first case we follow steps 9 
and 10 of Algorithm [21] to obtain the multiplication matrices of the radical with respect to its basis 
T=[l,x 1 ]: 
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with respective eigenvalues [2,-1] and [2,3]. 

For the second case we use the method described in [23] to compute the approximate multipli- 
cation matrices of the approximate radical of the clusters. Using Gaussian Elimination with complete 
pivoting, we found that the almost vanishing pivot elements were of the order of magnitude of 10 _1 
which clearly indicated the numerical rank. Using the submatrices obtained from the complete piv- 
oting algorithm we got the following approximate multiplication matrices of the approximate radical 
with respect to the basis T = [x\X2, £2]: 



The norm of the commutator of these matrices is 0.002 and their eigenvalues are respectively [1.949, —0.972] 
and [2.001, 3.036]. Note that the corresponding roots [1.949, 2.001] and [-0.972, 3.036] are within 10~ 2 
distance from the centers of gravity of the clusters, as was shown in [33J (recall that the radius of 
the clusters was 10 _1 ). 

In the third case, the numerical rank was not easy to determine using either SVD or complete 
pivoting. However, when we assume that the numerical rank of R is 2, and we cut the matrix R 
using the output of the complete pivoting algorithm, then we obtain the multiplication matrices with 
respect to the basis T = [x±X2, X2]' 



The norm of the commutator of these matrices is 0.010 and their eigenvalues are respectively [1.997, —0.987] 
and [1.999,2.993] (recall that the perturbation of the polynomials was of size 10~ 3 ). 

5 Conclusion 

In this paper we gave an algorithm to compute matrices of traces and the radical of an ideal X 
which has finitely many projective common roots, none of them at infinity and its factor algebra is 
Gorenstein. A follow-up paper will consider an extension of the above algorithm which also works in 
the non-Gorenstein case and for systems which have roots at infinity, as well as an alternative method 
using Bezout matrices for the affine complete intersection case to compute the radical \fX. 
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